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A new approach to global QCD analysis is developed. The main ingredients are two 
QCD-based evolution equations. The first one is the Balitsky-Kovchegov nonlinear equation, 
which sums higher twists while preserving unitarity. The second equation is linear and it 
is responsible for the correct short distance behavior of the theory, namely it includes the 
DGLAP kernel. Our approach allows extrapolation of the parton distributions to very high 
energies available at the LHC as well as very low photon virtualities, Q 2 <C 1 GeV 2 . 

All existing low x data on the F2 structure function is reproduced using one fitting 
parameter. The resulting x 2 /df = 1. 

Analyzing the parameter A = Slni^/c^ln 1/x) at very low x and Q 2 well below 1 GeV 2 
we find A ~ 0.08 — 0.1. A result which agrees with the "soft pomeron" intercept without 
involving soft physics. 
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1 Introduction 



In the paper JTJ a new approach to DIS was proposed. In the present paper we review and further 
develop the ideas introduced in Ref. 0. Our main result is that all existing low x data on the 
F2 structure function can be described by nonlinear QCD evolution. 

The standard perturbative QCD approach to deep inelastic scattering (DIS) is based on the 
DGLAP evolution equation || which provides the leading twist parton distributions. The main 
underlying assumption is that the high twist contributions are negligibly small if the evolution 
starts at sufficiently high photon virtualities ~ 2 — 4 (GeV 2 ). 

The approach based on the DGALP equation suffers from three principal problems. 

• The DGLAP evolution predicts a steep growth of parton distributions in the region of low x 
which would eventually contradict the unitarity constraints ||. Hence, we can expect large 
unitarity corrections to the DGLAP evolution equation, in the region of very low x. 

• The second problem is the general nature for any operator product expansion which is an 
asymptotic series. In application to DIS this means that the errors associated with the 
leading twist approximation are not small. They are of the order of the next to leading 
order twist contribution which grows very fast at low x. In fact, it can be shown that 
high twist contributions grow with decreasing x faster than the leading twist [|J . Hence, we 
cannot conclude that the higher twist contributions are small in the whole kinematic region, 
even if they are small for the initial value of Q 2 = Qq. The estimates of Ref. show that 
all available parameterizations of the solutions to the DGLAP evolution equation lead to 
substantial higher twist contributions. 

• The last problem is in the total nonability of the DGLAP equation to describe physics 
of low photon virtuality Q 2 < l(GeV 2 ). For these kinematics one needs to use Regge 
phenomenology or other phenomenological models. 

It is important to add that NLO DGLAP, though it improves the fits to presently available data, 
does not solve any of the above principal difficulties. Consequently, we are lead to the conclusion 
that DGLAP is insufficient to describe all kinematical phase space. For small values of x and/or 
Q 2 there is need for a new QCD-based idea. 

In this paper, we develop such an idea which allows us to extrapolate parton distributions to very 
low x (high energies). An extrapolation of the available parton distribution to the region of lower 
x is a practical problem for the LHC energies. We need to know the parton distribution both for 
estimates of the background of all interesting processes at the LHC, such as Higgs production, 
and for the calculation of the cross sections of the rare processes which are likely to be measured 
at the LHC. 

Our method was originally proposed in Ref. []J . It consists of two steps. As a first step, a nonlinear 
evolution equation, which takes into account the most significant higher twist contributions, is 
solved. This equation Q2.1| ) specifies a high energy (low x) behavior of the parton densities. The 
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solution obtained, below denoted as N, takes into account collective phenomena of high parton 
density QCD and respects unitarity constraints. Moreover, N can be found for large transverse 
distances, and so it also provides a possibility to describe data of low photon virtuality. 

The parton distributions which we obtain are then amended by adding to the solution of the 
nonlinear equation N, a correcting function AiV, which is aimed at correctly incorporating the 
short distance behavior of the theory, namely the DGLAP kernel. 

For AiV we propose a linear DGLAP-type evolution equation. The function AiV is considered to 
be a small correction to N concentrated in the region of moderate x. Consequently, this function 
should be free of all difficulties inherent in the usual solutions of the DGLAP equation. 

The philosophy of our approach is quite similar to the one recently presented by Ref. 0. That 
paper is a development of the Golec-Biernat - Wusthoff (GBW) model which in addition to the 
original model |7j] is improved by DGLAP evolution. We can trace a certain analogy between our 
function N and the original saturation model. Both functions play the very same role: they take 
into account the gluon saturation effects in unitarity preserving way, and describe physics of large 
distances. However, contrary to the saturation model the function N is derived from QCD. The 
DGLAP improvement both in our approach and in Ref. || is aimed at correctly incorporating 
the short distance dynamics however, the technical realizations are quite different. 

The paper is organized as follows. In the next section we review our approach and write the 
BK nonlinear equation for N and a linear equation for AiV. Section 3 presents some analytical 
estimates of the corrections induced by the DGLAP kernel. Section 4 is devoted to technical 
details of numerical solutions of the equations. The following Section (5) presents our fit to the 
experimental data on structure function and its logarithmic derivatives. Some predictions for 
THERA and LHC are given. Section 6 presents a general discussion which emphasizes some weak 
points of our approach. In the concluding Section (7) we summarize our results and mention our 
plans for future work. 



2 Reviewing a new approach to DIS 

The DGLAP equation describes the gluon radiation which leads to a strong increase in the 
number of partons. However, when the parton density becomes large, annihilation processes 
become important and they suppress the gluon radiation taming the rapid growth of the parton 
density || || ||. A development of new theoretical methods applicable to physics of high density 
QCD || ||, ||, [10], [11], [12], [L3|, |TJJ lead finally to the very same nonlinear evolution equation, 
nowadays credited to Balitsky and Kovchegov (BK)[|: 



Eq. (2.1) was originally proposed by Gribov, Levin and Ryskin p] in momentum space and proven in the 
double log approximation of perturbative QCD by Mueller and Qiu H. In the leading lnl/x appro xim ation it 
was derived by Balitsky in his Wilson Loop Operator Expansion (l2). In the form presented in Eq. (2.1) it was 
obtained by Kovchegov in the color dipole approach (l5| to high energy scattering in QCD. This equation was 
also obtained by summation of the BFKL pomeron fan diagrams by Braunfllq] and in the effective Lagrangian 
approach for high parton density QCD by Iancu, Leonidov and McLerran [|14| . Therefore, it provides a reliable 
technique for an extrapolation of the parton distributions to the region of low x. 
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2 l 2 iV(x 2, y; b - ^x i2 ) - iV(x 02 , b - *xi 2 )iV(xi2, y; b - ^x 02 " 



The equation is derived for N(r±,x;b) which stands for imaginary part of the amplitude of a 
dipole of size r± elastically scattered at the impact parameter b. 

In the equation (|2.1| ), the rapidity Y = — lnx and Y Q = — lnx . The ultraviolet cutoff p is needed 
to regularize the integral, but it does not appear in physical quantities. In the large N c limit 
(number of colors) Cf = N c /2 (we set N c = 3 in the numerical computations). 

Eq. ( |2.1| ) has a very simple meaning: the dipole of size x 10 decays in two dipoles of sizes x 12 and 
x 02 with the decay probability given by the wave function |\1/| 2 = 2 Xoi 2 . These two dipoles then 

X 02 X 12 

interact with the target. The non- linear term takes into account a simultaneous interaction of 
two produced dipoles with the target or, in other words, the Glauber corrections for dipole-target 
interaction. 



The linear part of Eq. ( |2.1|) is the LO BFKL equation ||17|| , which describes the evolution of the 
multiplicity of the fixed size color dipoles with respect to the energy Y. The nonlinear term 
corresponds to a dipole splitting into two dipoles and it sums the high twist contributions. Note, 
that the linear part of Eq. (|2.1|) (the BFKL equation) also has higher twist contributions and vice 
versa, the main contribution of the non-linear part is to the leading twist (see Ref. M for general 
arguments and Ref. for explicit calculations). 

As has been mentioned, the master equation (|2.1| ) is derived in the leading ln(l/x) approximation 
of perturbative QCD. This means we consider asin{l/x) « 1 while as <C 1 as well as ashiQ 2 <C 
1. In other words, the equation sums all contributions of the order (as ln(l/x))™ and neglects 
contributions of the orders as{as ln(l/x)) n and as In Q 2 (as ln(l/x)) n . Contributions of the latter 
will be taken into account by the function AN to be discussed below. 

Eq. (|2.1|) sums all diagrams of the order 



(a s (l/x) A ^J with A cc as ■ 



This means that starting from «5ln(l/x) ~ ln(l/«5) corrections due to rescattering and re- 
combination of parton become essential (see Ref. [0 for details). 

The next to leading corrections to the DGLAP or/and BFKL equations lead to A = Cias + C2a s 
which start to be important only for asln(i/a;) > 1/as- Therefore, the correct strategy is first 
to solve the master equation taking into account all corrections of the leading order, and only as 
a second step consider the next-to-leading order corrections. 
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It is well known starting with Bartel's paper in Ref. [Q that Eq. Q2.1|) can be proven in the 
large N c limit of QCD. Actually, this equation is a first theoretical realization of the Veniziano 
topological expansion For N c 1, we assume that a s = N c as/ir ~ 1 while aj C 1. 

An interesting feature of the equation is that it depends on as only, and all problems with the 
accuracy of the large N c expansion are really concentrated in the N c dependence of the initial 



distributions (see Ref. |18f for details). 

It should be stressed that a correct evolution equation without an additional assumption on large 
N c is known (see Weigert's paper in Ref. |ll|]). However, this equation is so complicated that we 
are far away from its solution. However, it is worthwhile mentioning that in the simplified case 
of double log approach (both asln(l/a;) and aglnQ 2 are of the order unity while as <C 1) the 
equation for N c m 1 was written and solved in Ref. ||20|| . This solution shows that corrections 
to large N c approximation are rather small and, therefore, at the first stage it is reasonable to 
neglect them. 



In order to safely use Eq. ( [2.1|) we neede to estimate the neglected contributions. A first class of 
such contributions is the interaction between two parton showers which leads to (as/N 2 ) ln(l/x) 
corrections which result in a bound on minimal x: 

N 2 

ln(l/z) < 

a s 

A second constraint comes from so called enhanced diagrams. It turns out that they lead to the 
very same restrictions as the previous one (see Ref. for details). 

The above energy limit is not very essential since the unitarity bound iV = 1 is reached at higher 
values of x. Thus the 1/N C corrections cannot modify this result, but could slightly modify the 
value of the saturation scale. 

The total dipole cross section is given by the integration over the impact parameter: 

<Xdipoie(r±,x) = 2 J d 2 b N(r ± ,x;b). (2.2) 

The contribution to the deep inelastic structure function F 2 which is due to N we denote by F 2 
and it is related to the dipole cross section 

F 2 (x,Q 2 ) = | d 2 r ± J dz P^(Q 2 ;r ± ,z) <x dipolc (r ± , x) . (2.3) 



The physical interpretation of Eq. (|2.3j ) is transparent. It describes the two stages of DIS [plf . 
The first stage is the decay of a virtual photon into a colorless dipole (qq -pair). The probability 
of this decays is given by P 7 *. The second stage is the interaction of the dipole with the target 
(c"dipoie in Eq. (|2.3D). This equation is a simple manifestation of the fact that color dipoles are the 
correct degrees of freedom in QCD at high energies fl5| . The QED wave functions of the virtual 



photon are well known |15, 22], g3[ (we consider only massless case) 



P^(Q 2 ;r ± ,z) 2 = ^ 2 J2Z 2 {(z 2 + (l-z) 2 )a 2 K 2 (ar ± ) + 4 Q 2 z 2 (1 - z) 2 K 2 (ar ± )} 



(2.4) 
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with a 2 = Q 2 z(l - z). 

It can be seen that Eq. ( [2.1Q does not depend explicitly on the target^. All the dependence on 
the target comes from the initial condition specified at some initial value Xq. For a target nucleus 
it was argued in Ref. [IS| that the initial conditions should be taken in the Glauber form: 



iV(x O i,x ;6) = N GM (x 01 ,x ;b) 



with 



N GM (x 01 ,x;b) = 1 - exp 



a S' 7rx 01 ^riDGLAP 



2N r R 2 



xG L 



(x,4/x^)5(b) 



(2.5) 



(2.6) 



The equation (|2.6|) represents the Glauber-Mueller (GM) formula which accounts for the multiple 
dipole-target interaction in the eikonal approximation \ 



24], [25]. The function S{b) is a dipole 



profile function inside the target. The value of xq is chosen within the interval 

/ 1 N 1 

exp( ) < x < — — , 

as 2mR 



(2.7) 



where R is the radius of the target. In this region the value of Xq is small enough to use the low x 
approximation, but the production of the gluons (color dipoles) is still suppressed as as ln(l/x) < 
1. Consequently, in this region we have the instantaneous exchange of the classical gluon fields. 
Hence, an incoming color dipole interacts separately with each nucleon in a nucleus (see Mueller 
and Kovchegov paper in Ref. [|I~0|1 ). 

For the hadron, however, there is no proof that Eq. ( [2.5|) is correct. Our criteria in this problem 
(at the moment) is the correct description of the experimental data. Almost all available HERA 
data can be described using Eq. ( |2.5| ) [p6| , |27j] , and we feel confident setting Eq. ( [2.5|) as an 
initial condition for Eq. ( |2.1|) . In our model, the Gaussian (S(b) = e~ b2 / R2 ) form for the profile 
function of the hadron is mostly used. The parameter R is a phenomenological input, while the 
gluon density xG DGLAP is a solution of the DGLAP equation. For a hadron target Eq. (|2.7| ) is 
still correct, but practically xq = 10~ 2 is chosen. This value satisfies Eq. ( |2.7| ) for which much 
experimental data exist, so one can check the initial conditions. 

Solutions to the BK equation were studied in asymptotic limits in Ref. P3 while several numerical 
solutions were reported in Refs. |16|, [1], |29], [50[ [31]]. In Ref. QTJ] and in present paper we solve Eq. 
( |2.1| ) in the coordinate representation in which the initial conditions are of a very simple form (see 
Eq. ( |2.5| )). The second reason for using the coordinate representation is the fact that all physical 
observables can be expressed in terms of the amplitude for the dipole-target interaction in the 
coordinate representation. Finally, it is also very useful that the long distance asymptotics is 
known: N — > 1 being N < 1 otherwise. This fact provides a natural control for the numerical 
procedure. 



Unfortunately, Eq. (24.) is an approximation. It only sums large In a; contributions. The situation 
can be improved at short distances. The exact x dependence of the kernel at short distances is 



2 This independence is a direct indication that the equation is correct for all targets (hadron and nuclei) in the 
regime of high parton density. 
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known, namely it is the DGLAP kernel. An attempt to obtain the elastic amplitude N based on 
elements of both the BK and DGLAP equations was presented in Ref. ||32j| . The authors of this 
paper first solve a generalized DGLAP-BFKL linear equation and then add to the solution 
a nonlinear perturbation of the form presented in Eq. ( f2.1| ). This approach actually incorporates 
the high twist contributions in the standard way, treating them as corrections to the leading one. 

We suggest a different approach to the problem. First, all twist contributions should be summed 
by solving Eq. (|2.1|) . Unfortunately, it is complicated to find a solution for an arbitrary value 
of the impact parameter b. We simplify the problem by solving Eq. ( |2.1| ) without including the 
6-dependence, this corresponds to case of solving for the initial condition at b = 0. At the very 
end we restore the 6-dependence by using an ansatz to be discussed in Section 4. We denote by 
N(r±,x) the solution of Eq. at b = 0. 

Secondly, we add to the solution obtained a correcting function AiV, which will account for the 
DGLAP kernel (Fig. [I]): 



N = N + AN for x < x 



[2.8) 



1/x 

A 



A N = 



evolution for A N 



N = N + AN 



BK 



N = N 



BK (evolution for N) 



A N = 



N 



DGLAP 



Figure 1: The kinematic map for 
the solutions. 



In order to extract the leading In (l/r^) contributions we define a set of new functions: 

n = N/{a s r 2 ± )- n = N/{a s r]_)- An = AN/{a s r 2 ± ). 

For the function n we propose the following nonlinear equation assumed to be valid in the leading 
large ln(l/r[) approximation: 

dn(r±,x) Cfo. s f 1 . . , x Cfo^t^ f 1 dz 

Pg^ g {z) n(r±, -) dz - 



d ln(l/rj 



7T Jx 



71 



J fn 2 (r ± ,f). (2.9) 

Jx/xq Z Z 



Here P g -+ g (z) stands for the usual gluon splitting function: 



W*) = 2 



+ 



+ z(l-z) + n±-%L) S {l-z) 



(2.10) 



6 



Note that we assume nonlinear effects to be of no importance for x > Xq. Eq. Q2.1|) can be 
rewritten in the large ln(l/r^) approximation as: 



dh(r±,x) dn(r±,xo) Cfoi s f 1 dz x. CVa 2 ?" 2 r 1 dz 



x . 



+ _^ (r±? _) _ _^( r± ,-). (2.11) 

a m(l/ rjj a ln(l/rfj tt J x z z n J x /x z z 

Subtracting Eq. ( |2.11| ) from Eq. (|2.9| ) and assuming AA" to be small compared to N, we derive 
the equation for An(r±,x): 



dAn(r ± ,x) C F a s f 1 . . , x 

a wm = L p ^ z) An{r ^ i } dz - (2 - 12) 

2 -^ } ^N(r ± , X -)An(r ± , X -) + ^ f (p^ 9 (z)- 2 -)n(r ± A)dz 

IT Jx/x Z Z Z 7T Jx/xo \ ZJ Z 

dn(r ± ,x ) C F a s x 

- a j h/ 2 v H / P 9^g{z) n{r ± ,-)dz. 

d ln(l/rj_) it Jx z 

Equation ( |2.12|) is a linear equation valid in the leading ln(l/r^) approximation. The first term on 
the right hand side of Eq. (|2.12|) is the DGLAP evolution for the correcting function AiV, while 



the second term is the "nonlinear interaction" of the solutions. The third term in the equation 
represents the correction which is due to the substitution of the BFKL kernel 1/z by the correct 
DGLAP kernel. 

The last term in ( [2.12j ) is the only term which accounts for the contribution of high x > Xq. In 
that region the function n = tt xG dglap / (2N C R 2 ) being a solution of the DGLAP equation is not 
given by a sum of n and An. The sign "-" in the upper integration limit indicates that in the 
limit x — > xq the 5-function term of the splitting function must be excluded. 

If we set A^ = then Eq. ( |2.12| ) reduces exactly to the gluonic part of the leading order DGLAP 
equation. In planned further development of our approach quark distributions and their evolution 
will also be included. At the present stage we take them into account implicitly in a s and by 
setting rif = 3 in P g -> g - 

The last two terms in Eq. ( |2.12j ) were omitted in Ref. |lj as they are not important at low x. 



However they should be included for correct computations. 

The initial condition AA^(rj_o,a;) = A^(rj_o,x) — A^(rj_o,x) is a phenomenological input at some 
initial transverse distance r_|_ to be specified. However, our strategy is based on the assumption 
that A^ describes the long distance physics correctly. Consequently, in order to eliminate any 
discontinuity of the function A^ we require 

AN(r ±0 ,x) = 0. (2.13) 

The value of r± is not specified but it is expected to be of the order 2 (GeV -1 ) corresponding to 
the naive relation Ql = 4/r^ ~ 1 (GeV 2 ). 
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Since we assume N(r±,x ) to describe correctly the data at x — x , the continuity across x = x 
supposes that iV(r_|_,x ) = N(r±,x ). If we require this equality then Eq. ( |2.12| ) would respect 
it provided initial condition (|2.13| ) is imposed. 

It is important to emphasize that Eq. ( ^.12j ) is free from the main problems of the DGLAP 
equation. First of all, the high twist contributions are summed (at least partially) by Eq. ( |2 . 1| ) . 
Secondly, our method respects unitarity. This is achieved due to the second term of Eq. ( f2.12p 
and the unitarity preserving initial condition ( [2.13| ). 



Finally it is necessary to compute a correction to F 2 structure function due to AiV. To achieve 
this goal we need to assign an impact parameter dependence of AiV. Similarly to what is common 
in DGLAP solutions, the 6-dependence is assumed to be a product of (2n R 2 )AN times a profile 
function. After ^-integration the latter contributes unity. Then the correction to F 2 reads 



2 



4 TV r r ± o c\r , 

AF 2 (x,Q 2 ) = —I / -± AN(t ± ,x)(ttR 2 ). (2.14) 

9 7T d JA/Q 2 r]_ 



The r_i_ integration in fl2.14|) is evaluated in the same large ml/r_i_ approximation as is valid for 



Eq. ( |2.12| ). Note that the coefficient (vri? 2 ) is cancelled with the same factor hidden in AiV. 



3 DGLAP correction - analytical estimates 

In this section we would like to make some comments regarding the consistency of our approach. 
It was argued previously that it is necessary to add a correction term AiV to the solution N of the 
nonlinear equation ( j2.1| ), where the function AiV is a solution of the evolution equation ( |2.12| ). 



Consistency of the approach requires the function AiV to give vanishing contributions to the 
dipole cross section at very small x. We also expect this function to decrease as r\ decreases. 
Finally, AiV is assumed to be a small correction to the function N. In order to check the above 
conditions some asymptotic estimates can be made without explicitly solving Eq. ( |2.12| ). Indeed, 
we will show below that Eq. (|2.12|) respects all the above mentioned requirements. 



Limit 1: fixed r±, x — > 0. 

At very small x and fixed distances the function N ~ 1. Eq. ( |2.12| ) can be simplified: 
am(l/r^J n Jx/x \ zJ \ z z J 



The main observation is that the evolution kernel entering the equation ( |3.15| ) is actually 
negative. Hence the function An decreases as r±_ decreases. 

Let us consider a model where the anomalous dimension has the form consistent with energy 



conservation 34 
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7 (w) =a s [--l), (3.16) 



where as = asN c /n and the anomalous dimension is defined by the Mellin transform of the 
splitting function P g ^ g - 



a s C, 



7 ( W ) = / dzP g ^ g {z)z". (3.17) 

7T JO 

Eq. ( p,15| ) can be solved using the inverse Mellin transform |TJ]. Define An(r_|_, cj) and 
u(rj_, a;) as the inverse Mellin transforms 

An(ri,w) = / dujx~ w An(r±, u) ; n(r_i_,u;) = / du n(r±,uj) . 

2m Jc 2ni Jc 

In the momentum representation Eq. (|3.15| ) together with the anomalous dimension ( |3.16|) 



is: 

= -*> ^ w ) + An ^ ^ ( 3 - 18 ) 



Eq. ( 3.18|) can be easily integrated. Applying again the approximation N ~ 1 we get the 



result for the correcting function AiV: 

AN(r x ,x) ~ — iV(r ± ,x) + C(x) (r 2 ± ) 1+ ° s . (3.19) 

l + a s 

The function C(x) should be determine from the initial condition AN(r±o,x) = 0. Conse- 
quently 

AN(r ± ,x) ~ - j^— \N(r ± ,x) - N(r ±0 , x) (rj_/r± ) 1+&s ] . (3.20) 

As expected the function AiV is negative and of the order 0(as) compared to N. As r± 
decreases, AiV decreases until it reaches a minimum at r± min determined by the equation 

^%T^\r,=r ±min - (l+a 5 ) N(r ±0 ,x)-^^,. (3.21) 
At shorter distances AiV tends to as it should. 



Note that at fixed r_|_, AiV is finite and non vanishing at x — > 0. Yet, this is consistent 
with the requirement of a vanishing contribution to the dipole cross section since the latter 
implies integration over the impact parameter b. We will assign different 6-dependences to 
the functions iV and AiV. After the integration, the dipole cross section due to iV will grow 
logarithmically with x, while the contribution of AiV will remain finite. 
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Limit 2: fixed x, r± — > 0. 

We would now like to address a question of the short distance asymptotics. In this limit, 
the function N is given by the solution of the BFKL equation. Namely, 

N(r±,u) ~ e° 5ln(1/ ^ )/w . (3.22) 
On the other hand, energy conservation (|3.16|) requires 



N(r±,u) ~ e ^sHi/r 2 ± ) (i/«-i) _ 

Hence 

AN(r ± ,uj) ~ iV(r ± ,u;) ( e MiAi) _ i) _> -N(r ± ,oj). (3.23) 

Finally we conclude that the function AiV is supposed to be negative. We expect AiV ~ (3 N 
with |/3 1 < 1 and to be approximately x independent. At short distances AiV tends to zero as N. 

Therefore, AiV turns out to be small in the whole kinematic region. The analysis presented above 
justifies the self consistency of our approach and paves the way for numerical calculations to be 
presented in the next section. 



4 Numerical solution of the equations 

In this section we report on the exact numerical solution of Eq. ( |2.1| ) with initial condition fl2.5|) 
and of Eq. (|2.12 ) with initial condition (|2.13 ). 



First of all we wish to discuss several technical details. 
• Kinematic domain 



The kinematic region where the solutions of fl2.1|) and ( 2.12|) are found, covers x values from 



x = Xq = 10~ 2 , where the initial conditions are set, to x — 10 -7 . The maximal transverse 
distance r_|_ is taken to be two fermi. The value of the ultraviolet cutoff p is 2 x 10~ 3 (GeV -1 ). 
The numerical solutions obtained are checked and are independent of this choice. 

• Coupling constant as- 

Eq. (|2.1|) is derived for constant as- However, the DGLAP equation and hence Eq. (|2.9|) have a 



running coupling constant. Consequently, the derivation of Eq. (|2.12|) implies the same running as 
in Eq. Q2.1|) . For numerical purposes we take the LO running as with as = as(4/r^) everywhere. 
At large distances we freeze as at the value as — 0.5. 

• Transverse hadron size R 2 . 

In Ref. ffl the fixed value R 2 = 10 (GeV -2 ) was taken. This choice corresponds to the value which 



is obtained from the "soft" high energy phenomenology |35|, |36|, and is in agreement with the 
HERA data on elastic J/^f photo-production (R7|. As R 2 is practically the only fitting parameter 



at our disposal we allow it to vary in order to fit the data. The optimal fit is achieved at the 
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value R 2 ~ 3.1 (GeV 2 ). This value is too small and requires our understanding. The physical 
meaning of such a small value will be considered in the Discussion. 

• Gluon density xG DGLAP . 

In our approach, the gluon density xG DGLAP appears twice: first, in initial condition (|2.5|) and, 
second, it accounts for the region x > x in Eq. ( f2.12j ). At this stage, we do not solve the DGLAP 
equation for the high x region. Instead, we rely on the existing parton distributions. Practically 
for xG DGLAP (x > Xo,4/r^) we use the LO CTEQ6 parametrization |5S . 

• Solution of Eq. (Q. 

In Ref. M Eq. ([2.1|) was solved by the method of iterations. In the present work we adopt 



another method which appears to be more efficient. Namely, we solve Eq. ( j2.1| ) as an evolution 
equation in rapidity with a fixed grid in r± space and a dynamical step in Y. The results of 
the new program are in total agreement with the old method of Ref. provided the very same 
initial input is used. The function N is shown in Fig. |2| (solid curves). At large distances, N 
saturates to unity, which is the unitarity bound. At short distances, N tends to zero indicating 
color transparency. 

x = 10 -6 





Figure 2: The solution of Eq. [2A_) (solid line) and Eq. ^2.1^J (dashed line) plotted versus r±. 



• Large distances. 

It is of crucial importance for our purposes to correctly determine N at large transverse distances. 



However, initial conditions ( |2.5| ) are not computable at large distances. The gluon parametrization 
appearing in Eq.(|2.5|) ends at r± ~ 0.5 fm. A resolution of the problem was suggested in Ref. 
39|. The function N possesses a property called geometrical scaling. Namely, N(r±, x) is not a 
function of two independent variables but rather a function of single variable t = r± Q s (x). Here 
Q s {x) stands for the saturation momentum scale. 



The solution N presented in Fig. ^] is obtained in two steps. First, initial conditions ( |2.5| ) are 
extrapolated to long distances by a constant, which at very long distances does not approach 
unity (such a procedure is not consistent with scaling that is a purely dynamical property of the 
evolution equation). Then a solution is obtained for all x. At sufficiently low x (x < 10~ 4 ) the 
initial conditions are forgotten and the dynamics is governed by pure evolution. In this region 
the geometrical scaling is manifested. As a second step, we take the solution thus obtained at 
x ~ 10~ 6 and use geometrical scaling in order to rescale this solution up to x = xq. The resulting 
curve is now used for a new long distance extrapolation of initial condition (|2.5|). The initial 
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condition obtained this way provides a smooth extrapolation of the Glauber formula to unity at 
very large distances. 

Finally, we note that the procedure presented above can be used for large distance extrapolation 
of the gluon density at high x > x$. 

• 6-dependence of the solution. 

We assume that solution of the equation Q2.1|) preserves the same 6-dependence as introduced by 
the initial conditions (12 - 



N{r x ,x;b) = (1 - e"^^), (4.24) 
where k is related to the 6 = solution 

K (x, r±) = - ln(l - N(r±, x,b = 0)). (4.25) 



A factorized form of the 6-dependence was recently advocated in Ref. ||40|| . The ansatz is quite 
good at moderate x, though it becomes worse at smaller x [|], |3(J. The overall uncertainty of the 
approximation can be roughly estimated not to exceed 10%-20%. 

We now proceed with the evaluation of the dipole cross section Eq. ( |2.2| ). Having assumed ( |4.24j ). 
the dipole cross section has the form 



Cdipoie = 2-kR 2 [ln(/c) + Ei(k) + 7] . 



(4.26) 



In equation ( |4.26| ) 7 denotes the Euler constant, while E\ is the exponential integral function. 
The expression (|4.26|) predicts the In k growth of the dipole cross section, which is in agreement 



with the conclusions presented in Ref. p8 



• Continuity at x = x$. 

In Section 2 we discussed the continuity of the function N at x = x . One of the ways for 
its realization is to require AN(r^,xo) = which is fulfilled when N(r±,xo) = N(r±,xo). 
However, N(r±, Xq) is given by GM formula ( |2.5| ) plus the large distance extrapolation while 
N ~ as r\ xG DGLAP . Formally they do not coincide though numerical differences are not signif- 



icant. Nevertheless, we decided to force the equality iV(r_i_,a;o) 
following changes were introduced in Eq. ( |2.12|) . 



N(r±,xo). To achieve this, the 



a) 



N 



GM; 



b) 



dn(r±,x ) 



C F a s 



71 



dz P g ^g(z)n( y r_ L ,x Q /z) . 



<91n(l/ri) 

The above changes are minor. Practically they affect only the long distance behavior of the theory 
for x > Xq, which is not significant. 

• Solution of Eq. ( |2?L2D . 

Having obtained the function iV we can search for the correction AiV. Eq. ( [2.12| ) is solved 
similarly to Eq. (|2.1|) , but with a fixed grid in rapidity and dynamical step in r±. The initial 
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conditions (|2.13|) are set at r_|_ = r± . The parameter r_|_ is an adjusting parameter to be 
determined from the optimal fit. The dashed curves in Fig. |2] show the correcting function AN 
corresponding to r± = 2 (GeV -1 ). 

The function AiV obtained displays all the qualitative properties deduced analytically. As ex- 
pected, starting from zero at rj_ the function AiV decreases until it reaches a minimum and then 
it increases to zero again at asymptotically short distances. The ratio \AN\/N increases per- 
manently during the evolution reaching about 70% at the edge of our kinematic domain. Below 
x ~ 10~ 3 this ratio is almost x- independent. 



5 Results 

5.1 Fitting strategy 

Low x F2 data is used to determine the parameters of our model. The experimental data for 
x < 10~ 2 is taken from ZEUS HI], 03], HI ||43|| , and E665 HJJ experiments. The overall number 



of points is about 345. We actually use the very same data as Ref. [gj. Statistical and systematic 
errors are added in quadrature. Whole data sets are allowed to be shifted within the overall 
normalization uncertainity. We use this freedom to shift the low Q 2 ZEUS data down by 2% and 
E665 data by 3%. The HI data were shifted up by 3%. 

Our fitting procedure is divided in two steps. First, recall that in our approach the function N is 
supposed to describe correctly the kinematical region of very low x and large r±. The only fitting 
parameter for N is R 2 . We vary R 2 in order to find an optimal fit for the F2 subset of data below 
Q 2 ~ 1 (GeV 2 ) (about 100 points). The resulting X 2 /df = 1.2 is achieved for R 2 = 3.1 (GeV~ 2 ). 

The function iV obtained by fitting low Q 2 subset of the data, is not capable of describing all data 
points. The fit to all points leads to x 2 /df > 3, which is not good. The reason for this mismatch 
is certainly due to the absence of the DGLAP kernel in the evolution of N. In order to solve this 
problem we switch on the DGLAP correction AiV which is our second step on the way to the 



optimal fit. To achieve this Eq. (|2.12 ) is solved. 



For AiV the only fitting parameter is the position r±o at which the initial conditions ( |2.13| ) are 



set. It appears, however, that variation of this position acts as a fine tuning parameter only. 
The optimal fit is realized at r_|_ = 2 (GeV -1 ) in total agreement with the underlying theoretical 
assumptions. 



5.2 Fit to F 2 data 

In this subsection we present results of the fit to the low x F2 data. The structure function F2 is 
given by a sum of three contributions: 

F 2 — F 2 + AF + F2 SQ , (5.27) 
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where the first two terms are given by Eqs. (|2.3| ) and ( p. 14 ). These terms take into account only 
the gluon contribution to F 2 . In fact, gluons are related to singlet quark distributions. The third 
term in (|5.27 ) takes into account contributions of non-singlet quark distributions: 



p NSQ _ \- i 



V 



(5.28) 



In the present stage of our research we borrow the valence quark distributions (q() from the LO 
CTEQ6 parametrization. It is important to note, however, that these distributions decrease with 
decreasing x and are of practically no significance below x ~ 10~ 3 . In future we plan to develop 
a fully self consistent approach without relying on any known parametrization. 

Our central results are presented in Fig. [5] for small Q 2 (a) and for large Q 2 (b). The solid line is 
the best fit obtained with resulting x 2 jdj = 1. The dashed line is a result obtained without the 
DGLAP correction AN (x 2 /df > 3). 
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Figure 3: Fit to the F 2 structure function. 



The quality of our fit is of the same level as of Ref. 
on a basis of a new scaling saturation model of Ref. |45 



A similar quality fit was also obtained 
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5.3 dF 2 /d{lnQ 2 ) 



The logarithmic derivative of F 2 with respect to In Q 2 is presented in Fig. |] at fixed Q 2 (a) and 
at fixed x (b). Only comparison with HI data is shown though similar ZEUS measurement 
exists as well ||46|| . Note that these experimental data were not take into account in the fitting 
procedure. 
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Figure 4: The logarithmic derivative Si^/SlnQ 2 . 



5.4 dlnF 2 /d{lnl/x) 

In this subsection we present our computation of A = 91nF 2 /(9(ln 1/x). A comparison with the 
HI data [|7| is shown in Fig. [5] at fixed Q 2 (a) and at fixed x (b). Fig. |G] presents our prediction 
for A at very low x and small values of Q 2 . At fixed Q 2 , A decreases with decreasing x tending 
to zero in agreement with the unitarity constrain. At Q 2 well below 1 GeV 2 and x ~ 10~ 6 , 
A ~ 0.08 ± 0.01. This value of A coincides with the "soft pomeron" intercept of the Donnachie 
and Landshoff model (DL) []35|| . It is important to stress that the result is obtained on a basis 
of perturbative QCD. The only nonperturbative input in our approach is the freezing of as at 
large distances. In fact, it was conjectured in Ref. ES| that the soft pomeron may appear in 
perturbative QCD due to freezing of as- 



15 



(a) 



(b) 



0.6 
0.4 
0.2 

0.6 
0.4 

'X 0.2 
t-i 0.8 
£ 0.6 
^ 0.4 
0.2 

^ 0.6 
0.4 
0.2 

0.6 
0.4 
0.2 





q a =90 



^4 




Q 2 =0.£ 



— i m i n i i i m i 

Q z =3.5 



Q a =1.5 



Q ! =15 



Q 8 =5 



Q s =2 



Q 2 20 



q 2 =6.5 



i 

Q 2 =2.5 



10-« 10- 3 10 a 10* 10 3 10 s 10-" io a io- ; 
X 



0.6 
0.4 
0.2 

0.6 

^ 0.4 

R 02 

T3 

\ 0.6 

£ 0.4 

T3 

0.2 

0.6 
0.4 
0.2 




x = 0. 00105 x=0. 00165 Jl x = 0.0026 





x-0. 000065 J_ x=0. 000105 1=0.000166 




20 40 60 20 40 60 20 40 60 

Q s (GeV 2 ) 



Figure 5: The logarithmic derivative A = dlni^/Sln 1/x. 




Figure 6: The logarithmic derivative 
A = d In F% /<9 In 1/x plotted at 
low Q 2 and very low x. 
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5.5 Prediction for F L at HERA 



In this subsection we present our prediction for the Fl structure function. The result is obtained 
on a basis of the function N (longitudinal part of F 2 ) (Fig. [7]). The function AN is obtained 
within the leading logarithmic approximation, and in this approximation it does not contribute 
to Fl- Note that at relatively high values of x ~ 10~ 2 — 10~ 3 our prediction can be slightly 
underestimated since the contribution of the valence quarks is neglected. 




5.6 Predictions for LHC and THERA 
5.6.1 Gluon density 

From the solutions obtained we can compute the gluon density. To this goal we rely on the 



Mueller's formula [22] which relates the density to the elastic dipole-target amplitude: 



^•^^/^Cf (5.29) 

Practically we integrate in x' up to xq and then add xG DGLAP \xq). Fig. [| presents a comparison 
between the gluon density obtained from Eq. ( |5.29| ) and xG CTE ® . At very low x a significant 
damping of the density can be observed compared to the DGLAP predictions. 
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5.6.2 F 2 



The obtained model allows extrapolation of the parton distributions to very high energies. Fig. 
|9] presents our predictions for THERA and LHC kinematics. 




6 Discussion 

6.1 The transverse hadron size R 2 . 

As was pointed out above the optimal fit is achieved for R 2 = 3.1 (GeV -2 ). Such a low value 
requires understanding and below we present several explanatory arguments. 
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1. First of all, the Glauber - Mueller formula of Eq. ( p.6|) can be used for a proton target 
with great reservations because contrary to the nuclear case large inelastic diffraction is 
present. This process not only has a considerable cross section but also a quite different 
impact parameter dependence corresponding to a small value of the radius. The effect 
of two different radii is in agreement with the HERA data on elastic and inelastic J/ty 
photo-production |49[| . In the simple additive quark model (AQM) there are two kinds of 
processes: rescattering of a dipole off one quark and rescattering due to interaction with 
two or even three constituent quarks. The admixture of the inelastic diffractive processes 
can be taken into account (see Ref. |)(J ) by effective decreasing of R 2 in Eq. (|2.6| ) from 
10(GeV~ 2 ) to 5(GeV" 2 ). 

2. Second, in Eq. ( |2.6| ) we used the Gaussian parametrization for 6-dependence. On the 
other hand, the data on J/ty production requires assuming the profile function of the form 

s{b) = Vr^ k ^ (6 ' 30) 

which corresponds to the power - like (dipole) form factor in momentum transfer represen- 
tation: ^ 

F dipole (t)= (l _ 1 _ RH)2 . (6.31) 

Fdi P oie{t) describe a system with the same radius R as the Gaussian form factor. In t- 
representation the latter looks as 

F exp (t) = e^ R2t . (6.32) 

Practically, the solution to the non-linear BK evolution equation is obtained at b = 0. Note 
that Sdi P oie(b = 0) = 2 Scaussianip = 0) and this difference can be interpreted as an effective 
decrease in the value of R 2 . In fact, a relatively good fit to the low x data can be obtained 
with the dipole profile function at R 2 ~ 4.5 (GeV~ 2 ). 

3. In Eq. ( |2.6| ) we use the following expression for the dipole - proton cross section: 

agjc 2 ^ 4 

Cdipoie = — 77 xG(x, ~y) , (6.33) 

N r r \ 



However, the expression Eq. ( |6.33|) is an approximation valid for small values of the anoma- 



lous dimension 7 only. The correct expression for the cross section was obtained in Ref. 



°an>oie{r^x) = j^n 2 J 0(x,/ 2 ) (1 - M!i ^ (6 . 34) 



with <p = dxG(x,l 2 )/dl 2 being an unintegrated gluon density The gluon density 

xG = xG DGLAP is a solution of the DGLAP equation 

xG{x,l 2 ) = — / dux- u ? (w)e l(u)hi2 . (6.35) 
2 m Jc 
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Substituting Eq. (|6.35|) into Eq. (|6.34j ) and performing the /-integration we obtain: 



0~ dipole 



4 as 7T 2 



N c Jc 2k i 



„2 \ !-7(w) 



1 - 7(cj) \ 4 



r(i + 7M) 
r(2 - 7 (w)) 



. (6.36) 



It turns out that Eq. ( |6.36| ) can be approximately rewritten in the very compact form: 



0~ dipole 



4 as vr z 



dr ± xG(x, —pr) . 



(6.37) 



Eq. ( |6.37| ) and Eq. (|6.33| ) are quite different (see Fig. (0)) which again can be taken into 
account by reducing the value of R 2 . 



dipole 




0.5 1 1.5 2 2.5 3 3.5 

r.CGeV- 1 ) 



Figure 10: The dipole cross section is 
plotted versus r± at x = 10 -2 
for Eq. ( 6.3%) (upper curve) 
and for Eq. ( \6.3% ) (I ower 
curve), as = 0.2 for this 
plot. 



6.2 Geometrical scaling and saturation scale 



We briefly discuss the issue of the geometrical scaling displayed by the function N. Namely 
N = N(t) with r = r±Q s (x). The phenomena of geometrical scaling for a solution of the BK 
equation was studied analytically in Ref 



2S 
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pqj and established numerically in Ref. 
Recall that in the extrapolation of the initial conditions to the very long distances we relied on 
the scaling property. The function N is displayed as a function of r in Fig. O. 

Three comments are in order. 



Fig. |TT| is obtained assuming Q s (x = 10 



l(GeV). 



Within 10% accuracy the scaling holds for r > 1. For smaller r there is a noticeable scaling 
violation depending on the value of x. In fact, more significant scaling violation is found 
in perturbative region compared to the results of Ref. 
be due to the difference between 



3£ 



This discrepancy is likely to 
a constant value of as = 0.25 was used 
It was argued in Refs. |53] that the 



in Ref. 

while the present work is done with a running 0:5 
running of a s provides an important source for scaling breakdown, when penetrating the 
region of perturbative QCD. 
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The x-dependence of the saturation scale can be investigated using the scaling property. 
Parameterizing Q s ~ x~ q we find q = 0.18 ± 0.02. This value is about half the size of 
previous estimates of Refs. j39fl . The latter were obtained with the constant as = 0.25 
and the decrease of q is doubtless due to running of a s . Indeed, in Ref. we showed that 
in the case of running as, the saturation scale grows much slower than the fixed constant 
case. It is certainly interesting to investigate the dependence of g on a s . 




Figure 11: Geometrical scaling. N ver- 
sus t = r± Q s (x). 



6.3 Comparison with the GBW model 

It was mentioned in the Introduction that the solution to the BK equation (N) in describing 
the low x data plays the same role as the GBW original saturation model. Consequently, it 
of interest to compare these two models. Fig. |T2] shows the dipole cross section Q4.26D plotted 
together with the one of the GBW model. Note that due to the impact parameter integration 
the dipole cross section ( f4.26|) grows with decreasing x (logarithmically) while the GBW model 
reaches a saturation value. 

As a function of r± the behavior of the curves in Fig. [12] is quite different. This is a numerical 
coincidence that after the r± integration these dipole cross sections (improved by the DGLAP 
corrections) lead to a good description of the very same data. 




r,(GeV-') 



Figure 12: The comparison between (Jdipoie ( 4-%°] ) (solid curve) and the dipole cross section of the GBW 
model (dashed curve). 
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6.4 Shortcomings of our approach 



We would like to list several shortcomings of our approach and indicate future steps for their 
elimination. 

• One of the theoretical difficulties is in the fact that AN — > const(x) < at low x and fixed 
b. In spite of the fact that this limit does not contradict the unitarity constraints, it looks 
very unnatural that the dipole amplitude does not reach the maxim possible value N+AN = 
1. Indeed, this fact is an artifact of our approximation, namely, of the oversimplified form 
of the non-linear term in Eq. (|2.9| ) in which 1/z should be also replaced by the full kernel 
P g ^ g (z). In future we plan to treat this problem. Here, we want to recall that after 
integration over b, J d 2 b AN becomes much smaller than J d 2 bN due to the logarithmical 
grows of the latter function of x. 

• Our results are based on the CTEQ parametrization, which enters our calculations through 
gluon distribution at x > 10~ 2 and valent quark distributions. 

We attempted to switch to another parametrization (GRV98 |55fl ) but failed to reproduce 
a very good fit (x 2 /df — 2.3). The main reason for this failure is because at Xq = 10 -2 
and very short distances the gluon of GRV is smaller than the CTEQ gluon, by about 10%. 
This difference cannot be practically eliminated by adjusting of our fitting parameter R. 
Yet, the treatment of the dipole cross section in the from presented by Eq. ( |6.37 ) is likely 
to improve the situation. 

At present, we have to conclude that our current results are parametrization dependent. 
This requires us to reconsider the problem by producing our own DGLAP fit of high x data 
which would also include the quark distributions. 

• One of the central uncertainities of our approach is in the impact parameter dependence of 
the function N. The ansatz (424]) is certainly not fully correct though it preserves the main 



properties of the 6-dependence. In our approach, the uncertainty due to this ansatz was 
partly hidden in the fitting of the effective target size R 2 . In order to eliminate this problem 
it is highly desirable to solve Eq. (|2.1|) including the full 6-dependence of the solution. 



7 Summary 

A new approach to DIS based on summation of high twist contributions in the leading lnx ap- 
proximation is developed. The first step implies solution of the Balitsky-Kovchegov nonlinear 
evolution equation. Secondly, a linear evolution equation for the correcting function which incor- 
porates the correct DGLAP kernel in the leading In Q 2 approximation is derived. It is important 
to stress that both the equations are based on QCD and derived in several approximations. 

The BK equation ( |2.1| ) is solved numerically by the method of evolution. The solution leads to a 
saturation of the function N at large distances. However, the dipole cross section obtained is not 
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saturated as a function of x. Due to the b integration it grows logarithmically with decreasing x 
in a contrast to the GBW saturation model. 



The DGLAP correcting function AN was found as a solution of Eq. Q2.12j ). In agreement with the 



analytical estimates this function contributes at moderate values of x and provides a correction 
to the main contribution due to N. 

As a main goal of this work, the low x F% data is fitted in the whole kinematic region both for 
small and large photon virtualities Q 2 . The resulting x 2 /df = 1. In order to achieve this result 
practically only one fitting parameter is used. This fit determined the optimal model which later 
is applied to compute logarithmic derivatives of F 2 . Several predictions for the THERA and LHC 
are presented. 

Analyzing A = dhiF2/d(hil/x) at very low x and small photon virtualities we found A ~ 0.08 
which coincides with the " soft pomeron" intercept of the DL model. It is important to stress that 
this soft pomeron occurs as an effective result of multiple hard (BFKL) pomeron rescattering. 
Except freezing of as, no soft physics is introduced in our approach. 

The obtained model opens a possibility to address many questions in high energy phenomenology. 
At present we are working on DIS off nuclei as well as more exclusive processes such as of J/ip 
production. 
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